% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Replication "Deconstructing the Yield Curve"
% Crump and Gospodinov (2024)
% Date: 26-JUL-2024
% Function for simulations based on VAR(1) specification
% Called by: simsVAR1.m
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function simOut = simsVAR1Fcn(s,data2,data,p,sp)

    rng(12345+s);

    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Simulate data
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    simfwd = simulateData(data2,data,p,sp,1);
       
    simyld = cumsum(simfwd,2)./cumsum(ones(size(simfwd)),2);
    simdret = zeros(size(simfwd));
    simdret(1:p.H,:) = NaN;
    simdret(p.H+1:end,2:end) = simfwd(1:end-p.H,2:end) - simfwd(p.H+1:end,1:end-1);
    simxret = cumsum(simdret,2);

    simX = simfwd(:,p.factorFwdMats);
    simX1 = [ones(sp.T,1),simX];

    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Calculate OLS estimator and t-stat
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    betaHatAll = NaN(numel(p.yldMats),size(data2.X1,2));
    betaSEHatAll = NaN(numel(p.yldMats),size(data2.X1,2));
    betaSEHatAll_LLSW_EWC = NaN(numel(p.yldMats),size(data2.X1,2));
    betaSEHatAll_LLSW_NW = NaN(numel(p.yldMats),size(data2.X1,2));

    for i = 1:numel(p.yldMats)
        tmp = nwest(simxret(p.H+1:end,i),simX1(1:end-p.H,:),p.H);
        betaHatAll(i,:) = tmp.beta;
        betaSEHatAll(i,:) = tmp.beta./tmp.tstat;
        tmp = nwest(simxret(p.H+1:end,i),simX1(1:end-p.H,:),sp.S_LLSW_NW);
        betaSEHatAll_LLSW_NW(i,:) = tmp.beta./tmp.tstat;
        tmp = ewc(simxret(p.H+1:end,i),simX1(1:end-p.H,:),sp.nu_LLSW_EWC);
        betaSEHatAll_LLSW_EWC(i,:) = tmp.betahat./tmp.that;
    end    

    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % BH Bootstrap
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [simfacload,simevals] = OrderedPCA(simyld, 'cov');
    Y1 = simyld*simfacload(:,1:3);
    simehat = simyld - Y1*simfacload(:,1:3)';
    sims2hat = trace(simehat'*simehat)/numel(simehat);

    bsOutBH = bootstrapVAR(Y1,[],sp.B,p.BH_B_BC,p.BH_B_BC,true,true);
    
    % Initialize Matrices
    bsBetaHatAll_CG = NaN(sp.B,numel(p.yldMats),size(data2.X1,2));
    bsBetaSEHatAll_CG = NaN(sp.B,numel(p.yldMats),size(data2.X1,2));
    bsBetaHatAll_BH = NaN(sp.B,numel(p.yldMats),size(data2.X1,2));
    bsBetaSEHatAll_BH = NaN(sp.B,numel(p.yldMats),size(data2.X1,2));
    bsBetaHatAll_Oracle = NaN(sp.B,numel(p.yldMats),size(data2.X1,2));
    bsBetaSEHatAll_Oracle = NaN(sp.B,numel(p.yldMats),size(data2.X1,2));
    bsBetaHatAll_FullOracle = NaN(sp.B,numel(p.yldMats),size(data2.X1,2));
    bsBetaSEHatAll_FullOracle = NaN(sp.B,numel(p.yldMats),size(data2.X1,2));

    for b = 1:sp.B
    
        % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % BH Bootstrap
        % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        bsYldBH = bsOutBH.bsY1(:,:,b)*simfacload(:,1:3)' + sqrt(sims2hat)*randn(sp.T,numel(p.yldMats));
        bsFwdBH = [bsYldBH(:,1) bsYldBH(:,2:end).*(ones(sp.T,1)*p.yldMats(2:end))/data.period-bsYldBH(:,1:end-1).*(ones(sp.T,1)*p.yldMats(1:end-1))/data.period];
    
        bsDretBH = zeros(size(bsYldBH));
        bsDretBH(1:p.H,:) = NaN;
        bsDretBH(p.H+1:end,2:end) = bsFwdBH(1:end-p.H,2:end) - bsFwdBH(p.H+1:end,1:end-1);
        bsXretBH = cumsum(bsDretBH,2);

        bsX_BH = bsFwdBH(:,p.factorFwdMats);
        bsX1_BH = [ones(sp.T,1),bsX_BH];
        
        for i = 1:numel(p.yldMats)
            tmp = nwest(bsXretBH(p.H+1:end,i),bsX1_BH(1:end-p.H,:),p.H);
            bsBetaHatAll_BH(b,i,:) = tmp.beta;
            bsBetaSEHatAll_BH(b,i,:) = tmp.beta./tmp.tstat;
        end
                
        % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Oracle Bootstrap
        % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        bsFwdFullOracle = simulateData(data2,data,p,sp,1);

        bsDretFullOracle = zeros(size(bsFwdFullOracle));
        bsDretFullOracle(1:p.H,:) = NaN;
        bsDretFullOracle(p.H+1:end,2:end) = bsFwdFullOracle(1:end-p.H,2:end) - bsFwdFullOracle(p.H+1:end,1:end-1);
        bsXretFullOracle = cumsum(bsDretFullOracle,2);

        bsX_FullOracle = bsFwdFullOracle(:,p.factorFwdMats);
        bsX1_FullOracle = [ones(sp.T,1),bsX_FullOracle];
        
        for i = 1:numel(p.yldMats)
            tmp = nwest(bsXretFullOracle(p.H+1:end,i),bsX1_FullOracle(1:end-p.H,:),p.H);
            bsBetaHatAll_FullOracle(b,i,:) = tmp.beta;
            bsBetaSEHatAll_FullOracle(b,i,:) = tmp.beta./tmp.tstat;
        end       

        % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % CG Bootstrap
        % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        simdretapprox = simfwd(1:end-1,2:end) - simfwd(2:end,1:end-1);
    
        varFwd = var1BC(simfwd(:,end),p.CG_B_BC,false);
        Zfwd = varFwd.VBC(2:end,:);
        Zinit = [Zfwd, simdretapprox];
        Z = [Zinit; Zinit(1:sp.M-1,:)];
        idxRand = randi(sp.T-sp.M,[ceil(sp.T/sp.M),1]);
        idxBootstrap = kron(idxRand,ones(sp.M,1)) + kron(ones(ceil(sp.T/sp.M),1),(0:sp.M-1)');
        idxBootstrap = idxBootstrap(1:sp.T-1);
        bsZ = Z(idxBootstrap,:);
        bsDret = [NaN(1,size(simfwd,2)); NaN(size(simfwd,1)-1,1) bsZ(:,2:numel(p.yldMats))];
        bsFwdN = filter(1 ,[1 -varFwd.PhiBC],...
               [0; varFwd.muBC+(bsZ(:, 1)-mean(bsZ(:, 1)))],simfwd(1,end));
        bsFwd = buildForwards(bsDret, simfwd(1,:), bsFwdN);
        %{
        % Check buildForwards
        bsFwd2 = NaN(sp.T,numel(p.yldMats));
        bsFwd2(:,end) = filter(1 ,[1 -varFwd.PhiBC],...
               [0; varFwd.muBC+(bsZ(:, 1)-mean(bsZ(:, 1)))],simfwd(1,end));
        bsDret2 = bsZ(:,2:numel(p.yldMats));
        bsFwd2(1,:) = simfwd(1,:);
        for n = numel(p.yldMats)-1:-1:1
            for t = 2:sp.T
                bsFwd2(t,n) = bsFwd2(t-1,n+1) - bsDret2(t-1,n);        
            end
        end
        %}
        
        bsFwdCG = bsFwd;
        bsDretCG = zeros(size(bsFwdCG));
        bsDretCG(1:p.H,:) = NaN;
        bsDretCG(p.H+1:end,2:end) = bsFwdCG(1:end-p.H,2:end) - bsFwdCG(p.H+1:end,1:end-1);
        bsXretCG = cumsum(bsDretCG,2);

        bsX_CG = bsFwdCG(:,p.factorFwdMats);
        bsX1_CG = [ones(sp.T,1),bsX_CG];                

        for i = 1:numel(p.yldMats)
            tmp = nwest(bsXretCG(p.H+1:end,i),bsX1_CG(1:end-p.H,:),p.H);
            bsBetaHatAll_CG(b,i,:) = tmp.beta;
            bsBetaSEHatAll_CG(b,i,:) = tmp.beta./tmp.tstat;
        end
    end

simOut.betaHatAll = betaHatAll;
simOut.betaSEHatAll = betaSEHatAll;
simOut.bsBetaHatAll_CG = bsBetaHatAll_CG;
simOut.bsBetaSEHatAll_CG = bsBetaSEHatAll_CG;
simOut.bsBetaHatAll_BH = bsBetaHatAll_BH;
simOut.bsBetaSEHatAll_BH = bsBetaSEHatAll_BH;
simOut.bsBetaHatAll_FullOracle = bsBetaHatAll_FullOracle;
simOut.bsBetaSEHatAll_FullOracle = bsBetaSEHatAll_FullOracle;
simOut.betaSEHatAll_LLSW_EWC = betaSEHatAll_LLSW_EWC;
simOut.betaSEHatAll_LLSW_NW = betaSEHatAll_LLSW_NW;
simOut.simEigvals = simevals;
